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Graphene is a recently discovered carbon based material with unique physical properties. This 
is a monolayer of graphite, and the two-dimensional electrons and holes in it are described by the 
effective Dirac equation with a vanishing effective mass. As a consequence, electromagnetic response 
of graphene is predicted to be strongly non-linear. We develop a quasi-classical kinetic theory of 
the non-linear electromagnetic response of graphene, taking into account the self-consistent-field 
effects. Response of the system to both harmonic and pulse excitation is considered. The frequency 
multiplication effect, resulting from the non-linearity of the electromagnetic response, is studied 
under realistic experimental conditions. The frequency up-conversion efficiency is analysed as a 
function of the applied electric field and parameters of the samples. Possible applications of graphene 
in terahertz electronics are discussed. 
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I. INTRODUCTION 

The region of the electromagnetic spectrum from 0.3 to 20 terahertz (so called Terahertz gap) is a frontier area 
for research in physics, chemistry, biology, material science and medicine^. Due to the recent progress in THz 
technologjsii^, the THz studies continue to expand, involving more and more scientists in the development of new 
sources and detectors of THz radiation, as well as in research on applications of THz waves in different areas. The 
search for new methods of THz emission and detection and development of simple, compact and inexpensive THz 
sources and detectors remains a challenging problem. 

The most common technique for producing low-power electromagnetic radiation at frequencies above 0.3 THz is 
through nonlinear multiplication (upconversion) of lower frequency oscillators^, fl — *■ mfi, m ~ 2, 3, . . .. Such 
upconverters, which are commonly based on GaAs Schottky barrier diodes, successfully work as doublers and triplers 
(m — 2 and 3), but provide very poor conversion efficiency for higher order harmonics (m > 3)^'-. The search for 
alternative nonlinear materials, which could provide efficient frequency multiplication, especially with the upcoversion 
factor TO > 3, is therefore highly desirable. 

In this paper we discuss the recently predicted^ effect of the frequency multiplication in graphene. Graphene is 
a new material, which has been experimentally obtained about three years ago^^!^'^ and has immediately attracted 
great interest of researchers, for recent reviews see^i^^. Graphene is a monolayer of graphite, in which carbon atoms 
are packed in a dense two-dimensional (2D) honeycomb lattic o^-'^'^^i^^ . This is a two-dimensional semimetal, with 
a very specific electronic band structure. Fig. [TJ The lower {I = 1) and upper {I — 2) branches of the energy 
spectrum Epi touch each other at six points at the corners of the hexagon-shaped Brillouin zone. In the ideal case of 
a uniform, undoped graphene at zero temperature, the lower band Epi is fully occupied while the upper band Ep2 is 
empty, and the Fermi level goes through these six, so called Dirac points P*^*^ = TiK^'^ z = 1, . . . , 6. Near the Dirac 
points, at p ~ ?iK^'^\ the electrons in graphene have a linear, quasirelativistic dispersion with zero effective mass of 
quasiparticles, 

= {-iyv\p - hK(^^\ = {-lyv^pi+pi, (1) 

where p = p — h'K^'^\ Only two of the six Dirac points in the Brillouin zone, usually referred to as K and K' points, 
are inequi valen t 1 1 1 1 ^ 1 1 3 ^ ^j^g phenomena, discussed here, the states of the K and K' cones give additive contribution 
to the ac graphene response, so that one can consider only the states in the vicinity of one Dirac cone, accounting for 
the second cone by introducing the valley degeneracy factor — 2. The velocity y in ([T|) is about 10^ cm/s, according 



FIG. 1: The band structure of electrons in graphene (side view and top view). The two branches of the energy spectrum Epi, 
I = 1,2, touch each other at six points at corners of the hexagonal Brillouin zone. The values on the axis are k^a and kya, 
where a is the lattice constant of graphene. 
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to measurementsi!^. Near the K^*^ points the graphene quasiparticles are described by two-component spinor wave 
functions, determined by the effective Dirac equation, and are called Dirac fermions. 

The unusual, linear dispersion of graphene quasiparticles near the Fermi level leads to a number 
of interesting and so far not fully understood transport phenomena, such as the minimal electrical 
conductivity^'^ i-^'^'-^^i-^^'-^'^i-^^'-^^'^^i^^'^^i^^'^'^'^^ , absense of the weak localization^^, unconventional quantum Hall 
effecl!^i2^i21i22ii2S, observable even at room temperature^^, and other. Electrodynamic properties of graphene, which 
have been studied both experimentally^^^ and theoreticallj^^>i^^^>29^^^^^^^i^^i^^^ia, also 
demonstrate non-trivial features in the frequency dependent conductivityi^i2ii22ii2^ii^, photon-assisted transport^^, 
microwave and far-infrared respons o^^'^^i'^^'^^i"^°'^^ , plasmons spectrum^i^i^i^i^, et cetera. New electromagnetic 
modes, specific only for the graphene system, have been also predicted^i^. 

Both the transport and electrodynamic properties, briefly outlined above, have been studied within the linear 
response theory. Going beyond the linear-response approach, one can show^*^ that graphene should demonstrate 
strongly nonlinear electromagnetic response at relatively low amplitudes of the external electric field. In particular, 
irradiation of graphene by an electromagnetic wave with the frequency f2 should lead to the frequency multiplication 
f2 — > mf2 with odd values of to = 3, 5, 7, . . .. This makes graphene a simple and natural frequency multiplier— and 
opens up exciting opportunities for using graphene in terahertz electronics. 

In general, prospectives for building graphene based devices for terahertz applications are very attractive. Apart 
from the frequency multiplication, graphene devices could be used in p/asmon-based voltage-controlled sources and 
detectors of THz radiation. The physics of such devices has long been discussed in the literatur o^^'^-'^i^^i^'^i^'*'^^ and 
the main ideas of 2D plasmon-assisted detection and generation of radiation have been confirmed in a number of 
experiments on conventional 2D electron system o^^'^'^i^^'^^ . Graphene has that advantage that the Fermi velocity of 
electrons in it is much higher than in other semiconductor materials, and its plasma frequency is widely tunable and 
lies in the THz rang o^^i^^ . It should be noticed that a closely related to graphene material - carbon nanotubes - has 
also demonstrated a great potential for terahertz application a^°i^^'^^i^^'^'*i^^ . 

In this paper we discuss the effects of the non-linear electromagnetic response and the frequency multiplication in 
graphene^i^. We consider response of graphene to a strong uniform external time-dependent electric field within the 
semi-classical kinetic theory (Section |n|. In addition to results obtained in Refsi^i^, we take into account the sclf- 
consistent-field effects (Section III Cp . which should be important under realistic experimental conditions, and discuss 
both harmonic (Section |TT|) and pulse excitation (Section IIII|) of the system. The results obtained are summarized 
and discussed in Section HVl 



II. SELF-CONSISTENT KINETIC THEORY OF NON-LINEAR ELECTROMAGNETIC RESPONSE OF 

GRAPHENE 



A. Qualitative consideration 

Due to the linear dispersion ([T]) response of graphene to an external electromagnetic field turns out to be intrinsically 
non-linear, which naturally leads to the frequency multiplication effect. Physically, the possibility of the frequency 
upconversion in graphene can be explained in a very simple manner. Consider a classical 2D particle with the charge 
— e and the energy spectrum 

i?p2 = Vp = V^pI+pI (2) 

in the external time-dependent harmonic electric field Ex{t) ~ Eq cos fit, Figure[2^. From now on wc will consider only 
electrons in the vicinity of one Dirac cone, omit tilde in ([1]) and take into account the presence of two inequivalent 
cones in the Brillouin zone by the valley-degeneracy factor — 2. According to the Newton equation of motion 
dpx/dt = —eEx{t) the momentum Px{t) is then given by the sine function 

Px{t)=po{t)^~{eEQ/^l)sixint, (3) 

Figure [213. In conventional 2D electron systems with the parabolic energy dispersion, the velocity Vx, and hence, the 
current jx = —eUgV^ are proportional to Px, so that the normal 2D system responds at the same frequency O (here 
Us is the areal density of particles). In graphene, however, the velocity 



(4) 
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is a strongly non-linear function of px, therefore response of graphene is substantially anharmonic, Figure [2fc- In the 
extreme limit, when py in Eq. ([¥]) is close to zero, is proportional to sgn(pa;) and the ac electric current is 

jx(t) = ensVsgn(smQt) — eUgV— jsinilt + - sinSili + - sinSOt + . . .1 . (5) 

TT |_ 3 5 J 

The current ^ contains all odd Fourier harmonics, with the amplitudes jm, = 1, 3, 5 . . ., falling down very slowly 
with the harmonics number, jm ~ l/l'^^l- An isolated graphene layer should thus work as a simple and natural 
frequency multiplier, with the operating frequency, variable in a broad range. 



B. Kinetic approach 

The current ([5]) above is independent of the electric field, which means that Eq. ([5]) is not completely correct. The 
reason is that in the above simple consideration we did not take into account the Fermi distribution of electrons over 
the quantum states in graphene. To do this, we useS. the kinetic Boltzmann theory, which allows one to get an exact 
response of the system not imposing any restrictions on the amplitude of the external electric field E'^^*(i). 

In a real experimental situation, graphene sheet lies on top of a silicon oxide - silicon structure, and the gate voltage 
Vq can be applied between graphene and the silicon substrate, in order to control the density of electrons or holes 
in graphene. In addition, the graphene-Si02-Si system can be doped by impurities. Both the gate voltage and the 
doping can shift the chemical potential /i of electrons in graphene to the upper E-p2 or to the lower E^i band. Assume 
that the chemical potential /i lies in the upper band E^2 = the temperature is small, T <C and the system 
is subjected to the external time-dependent ac electric field E'^^*(<). Then the momentum distribution function of 
electrons /p(i) is described by the Boltzmann equation 



dt 



dp 



eE'="=*(i) = 0, 



(6) 



in which we have ignored collisions of electrons with impurities, phonons and other lattice imperfections. Equation 
ISl) has the exact solution 



/p(i) =-^o(p-PoW) 



(7) 



where 



^o(p) = 



1 + exp 



Vp- 

T 



(8) 
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is the Fermi-Dirac function, and 

Po{t) = -e f ^'^\t')dt' (9) 

J — OO 

is the solution of the single particle classical equation of motion. The electric current 

j(t) = -^E-/pW (10) 



s 

p 



then assumes the form 



9sgveV f pdp 

where gr^ = 2 is the spin degeneracy and S is the sample area. If the temperature is zero, T — 0, and the chemical 
potential is finite, > 0, the current j{t) can be rewritten in the form 

m = en^V ^_ g(Q), (12) 



where P = P(t) = —po{t)/pp, P{t) = |P(i)|, Pf — f^/V is the Fermi momentum, and 



(13) 



I Li a I Lp 1-1 r\ 

Anh^ 47rnV2 

is the density of electrons in the upper band. The function G{Q) in (fT2|) is defined as 

4 2P(t) 

GiQ)^^ cos xdxiy'l + QCOSX- cos x), Q{t) = ^' < 1, (14) 

and analyzed in Appendix 1X1 

If the external field E*^^' is small, so that the Fermi distribution is weakly disturbed, P{t) ^ 1, the function 
G{Q) « 1, see Eq. and 

j(t) ^ IlL^ / ^--\t')dt'. (15) 

VF 7_oo 

In this, linear response regime, the current is directly proportional to the electric field. If the external excitation is 
harmonic, 

E^^*(t) =EoCOsm, (16) 

the current is 

j(i)«^^Eosinm, (17) 



which corresponds to the linear-response frequency-dependent coUisionless scalar Drude conductivity 



24.28.29.37.38.48 



a{Vl) = — = 1—— — — . (18) 

vipF h An hn 

If the external field E*^^' is large and the Fermi distribution is strongly disturbed, P{t) can be much larger than 
1, the function G{Q) ~ 1, and response of the system to a harmonic excitation (|16p is described by Eq. ([5]). The 
non-linear regime ([5|) is thus achieved at \po \ ^ Pf, or at 

the energy, gained by electrons from the external field during the oscillation period should be large as compared to 
their average equilibrium (Fermi) energy. 
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FIG. 3: (Color online) (a) The time dependence of the ac electric current, measured in units enaV, at harmonic excitation of 
the system at the frequency Q. The temperature is zero, T//i = 0; the curves are labeled by the values of the electric field 
parameter £ = eEoV/Qn. (b) The Fourier components of the current (I12p as a function of £ at T/fi = 0. At f t ^ oo the 
curves tend to the values 4/nm, see Eq. ([5]). 



Figure [3^ illustrates the time dependence of the current (|12p at the harmonic excitation of the system (fTB]). One 
sees that in the low- field limit the response is linear (the j{t) dependence has a sinusoidal shape), while at strong fields 
the time dependence of the current tends to that given by Eq. ([5|) . In Figure [8)3 wc show the Fourier components 
of the current versus the field parameter £. One sees that when £ becomes larger than ~ 4, the Fourier amplitudes 
saturate and one gets in the ultimate non-linear regime. 

The strong-field condition (fT9|) can be rewritten as 

» (20) 

which shows that the required ac electric field grows linearly with the electromagnetic wave frequency and with the 
square root of the electron density. At / ~ 100 GHz and — 10^"'^ cm~^, the inequality ((20)l is fulfilled at > 200 
V/cm. 

The results obtained above refer to the case T / ^ = {) {ji'is finite). The opposite limit /i = (at a finite temperature, 
fi/T = 0) is difficult to realize in the whole sample, at least, in currently available systems. Due to inhomogeneity 
of the samples, the zero-energy point fluctuates in space, with a typical local electron or hole density of order of 
10^^ cm~^, Refi^. The condition ji/T = can thus be satisfied only locally, at certain points inside the sample. 
Nevertheless, for the sake of completeness (and in view of possible availability of higher quality samples in future), we 
also present here results^ for the case ^/T = 0. For arbitrary values of /z/T some results can be also found in Ref."'-. 

At a finite temperature and the vanishing chemical potential fi — both electrons and holes contribute to the 
charge carrier density 

Us = Tie + nil ~ ^ (21 

and to the current. Starting again from Eq. pip but accounting for the hole contribution and setting ^ = 0, we get 

j{t)=ensV^^ f xdx ( de . , "^""^^ (22) 

Jo Jo 1+ exp (y/x^ + - 2xPtCOs9] 

where the momentum po (t) is now normalized by the characteristic temperature-dependent momentum, Pt = Pt (t) = 
-Po(i)/PT, PT{t) = |PT(i)l: and pt = T jV . 

If the external field is small, Pt(J-) ^ 1, Eq. (|22p gives the linear-response current 

,,61n2 7ise2y2 61n2 /"* , , , 

j(t) « ensWrit)^ = ^ / E^-*(t')rft', (23) 
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FIG. 4: (Color online) (a) The time dependence of the ac electric current, measured in units eUsV , at harmonic excitation of 
the system at the frequency f2. The chemical potential is zero and the temperature is finite, ji/T = 0; the curves are labeled 
by the values of the electric field parameter £t = eEoV/QT. (b) The Fourier components of the current (|22|) as a function of 
the field parameter £t = eEoV/flT at fj,/T = 0. At f t — > oo the curves tend to the values A/nm, see Eq. (O. 



proportional to the electric field. If the external excitation is harmonic, Eq. (jl6p . the current is 

i(t) = ^J^ ^Eosinm, (24) 

which corresponds to the linear-response frequency-dependent coUisionless Drude conductivitj«2i 

61n2m,e2y2 M2e^g,g,T 

In the strong field regime P^o <^ 1 Eq. is reduced, again, to The dimensionless field parameter in the case 
/i/T = is defined as £t = eEoV/TQ. 

Figured^ demonstrates the time dependence of the current ([^^ at the harmonic excitation of the system The 
j{t) dependence is qualitatively similar to that, shown in Figure [3^. In Figure we show the Fourier components of 
the ac electric current, for to = 1, 3 and 5, as a function of the field parameter £t bA, fi = 0. 

In the case /^/T = the strong- field condition £ ^ I assumes the form 

Eo > (26) 
eV 

At the temperature T ~ 100 K and the frequency / ~ 100 GHz this corresponds to the requirement Eq > 50 V/cm. 

Using the quasi-classical (Boltzmann) approach for the description of electromagnetic response of graphene imposes 
certain restrictions on the validity of the presented theory. Physically, using the Boltzmann theory, one takes into 
account the intra-band contribution to the ac electric current, but ignores the inter-band contribution due to transitions 
between the lower (quasi-hole) and the upper (quasi-electron) bands. This is possible if the frequency of the external 
radiation satisfies the inequality 

hn ■^max{^l,T}. (27) 

In practically interesting cases (the density ~ 10^^ — 10^^ cm^^, room temperature) this restricts the frequency of 
radiation by ~ 10 — 30 THz. The presented quasi-classical theory can thus be used in the whole Terahertz gap. 



C. Self-consistent- field effects and radiative decay 



In Section IIIBI we have not considered effects of the radiative decay, which can be important under realistic 
experimental conditions in graphene. It was assumed that the graphene electrons move in the sample under the 
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action of the external electric field E^^*(i), and this directly leads to the time-dependent electric current Eq. 
P^ . In general, however, the time-dependent electric current creates, in its turn, a secondary (induced) electric 
field E""'(i), which acts back on the electrons and should be added to the external field. Calculating response of the 
system, one should take into account that electrons respond not to the external, but to the total self-consistent electric 
field E'°*(i) — E'^^*(t) -I- E*"'^(i). This results in the effects of electromagnetic reaction of the medium (graphene) to 
the external field and can reduce the frequency upconversion efficiency. How substantially suppresses the radiative 
decay the efficiency of the frequency multiplication, is studied below. 

Consider an infinite 2D electron system with the graphene sheet lying at the plane z — 0. We assume that the 
external electromagnetic wave is incident upon the graphene layer along the z axis and induces the ac current in the 
layer. This current produces the induced electric field E™'^(t), which is added to the external one. The Boltzmann 
equation for the momentum distribution function of the electrons should then be written as 

5/p(t) duit).^, 



dt dp 



eEltoit) = 0, (28) 



instead of ([5]). The solution of this equation, as well as the electric current, can again be written in the form ([7]) and 
(jlip . respectively, but the classical momentum po(t) now satisfies the equation 

Po(i) = -e /* El%it')dt' = -e [K%it') + K%it')] dt' , (29) 

where the field E*°*q(/:) is not known and should be calculated self-consistently. To do this, we recall that the current 
and the induced electric field are related by the Maxwell equations, 

ErioW = -27rj(t)/c. (30) 

Combining (P^ . ([50]) and pT|) we get the following self-consistent equation of motion for the momentum po(i): 

dpo(t) ^ e^gsg^V f pdp , 



dt 2TTh'^c 



-^o(p-Po(t)) = -eE-*o(i). (31) 



After the non-linear equation ([3T|) is resolved with respect to the momentum po{t), the current j(t) can be found from 
Eq. dm): 

gsQveV f pdp 

JW - /o-t;N? / —^0 (P - Po(0) : (32) 



Equations ((3T|) and (|32p describe the non-linear self-consistent response of graphene to an arbitrary external time- 
dependent electric field E|^*Q(i). The second term in the left-hand side of Eq. (|31|) describes the radiative decay 
effects in the infinite 2D graphene layer. 

In conventional 2D electron systems with the parabolic energy dispersion and the effective mass m* of 2D electrons 
the self-consistent equations for po(i) and j(t), similar to ([3T|) and ((32|) . have the form 

^ + ^PoW = -eErJo(t), jW=-^Po(t). (33) 
at m*c ni* 

Here Tpar = 27rnse^/m*c is the radiative decay rate^ in the conventional (parabolic) 2D electron systems. In ([55)1 
we have ignored the scattering due to impurities and phonons [the corresponding term 7Po(i) can be added to the 
left-hand-side of the first Eq. (l33|) ]. In high-electron mobility GaAs/AlGaAs quantum- well samples the radiative 
decay Tpar substantially exceeds the scattering rate 7, Tpar ^ 7, and determines the linewidth of the cyclotron, 
plasmon, and magnetoplasmon resonances^^iS^. As the graphene mobility is also very high, one can expect that at 
high frequencies the radiative effects are more important in graphene than the scattering effects. This justifies ignoring 
the scattering terms in Eq. ([5T|) . 

Returning back to the non-parabolic graphene system, we rewrite pip (at T = 0) in terms of the dimensionless 
momentum P(t) = —po{t)/pF- 

where 

g^gv 2fi 
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FIG. 5: The time dependence of (a) the momentum P{t) and (b) the current j{t)/ensV at the field parameter f = 10 at several 
values of the radiative decay T/Q. 



The current is determined, again, by Eq. 

In the linear-response regime, when \P\ <C 1 and « 1, Eq. (p4)) gives 

^ + rPW = -E-yi). (36) 

at pf 

From here one sees that the quantity T has the physical meaning of the radiative decay rate in graphene in the 
linear-response regime. In contrast to Tpar, T is proportional to the square-root of the charge carrier density. For 
experimentally relevant densities Ug the value of T lies in the subtcrahertz range, r/27r(TIIz)ss 0.13-\/ns(10^^/cm-^). 

Now consider response of graphene to a harmonic excitation ([TB|). It is convenient to rewrite Eqs. ([M)) and in 
the form 

dP r P , , eEnV , , 

from which one sees that the solutions depend on two dimensionless parameters, £ = eEoV/ jiVl and F/il. If the field 
parameter £ is small, the response is linear, \P\ <^ 1, and 

P(t)^— ^^cos(r^.-arctan^), (39) 

i{t)^en,V¥{t). (40) 

More interesting is the case of the strong electric fields, £ ^ 1. Figure [5] shows the time-dependence of the momentum 
P{t) and of the dimensionless electric current i{t)/ eUgV at £ = 10 and several values of the radiative decay parameter 
T /VI. If F/57 does not exceed the value of about £ /2 (=5 in the example of Fig. O, the self-consistent-field effects 
lead only to the phase shift of the current, not influencing the shape of the current-time curves and hence not reducing 
the amplitudes of the higher harmonics. At higher values of F/fi (between F/fi ^ £ /2 and F/f2 ~ £) the shape of 
the current-time curves smoothly modifies from the step-like form to the sinusoidal form, and at F/f2 ^ £ the higher 
harmonics are fully suppressed. Quantitatively this can be seen from analysis of the Fourier harmonics of the current. 
Expanding the current in the Fourier series 

- [A,n cos(mm) -f Bm s\n{mVlt)] , (41) 



9 



m= 


=1 


n- 


=3 


m= 


=5 


m= 


=1 


m= 


=9 




FIG. 6: Fourier harmonics of the current j{t)/ensV at f = 10 as a function of T/Q,. 



we calculate the amplitudes F„i ~ y/ A^^ + and plot them as a function of T/ft in Figure [6] for the lowest 
harmonics numbers from 1 to 9. One sees that, indeed, the higher harmonics are suppressed at F/fi > £, but almost 
not influenced by the radiative decay if F/f2 < 0.75. We have performed similar calculations at other values of the 
field parameter £ and found the same results if 5 ^ 1. 

The condition of the effective frequency upconversion, with account of the self-consistent-field effects and the 
radiative decay, thus assumes the form £ > F/il, or 



fiT 2TTnseV 
eV ~ c 



300— X n,(10"/cm2). 
cm 



(42) 



This condition does not depend on the frequency of radiation. 

It should be noticed that in the above consideration the sample was assumed to be infinite. Quantitatively, this 
means that the lateral dimensions of the graphene layer L should exceed the wavelength of radiation A. At the 
frequency of 1 THz the above derived formulas will thus be valid for samples with the dimensions larger than ~ 300 
/im. As the currently available samples are smaller (~ 10 — 30 /Ltm), the above consideration overestimates the role of 
the radiative decay. If L ^ A, the radiative decay rate F in the above formulas should be replaced by F(L/A)2, Ref.— . 
The upconversion effect should then be observed at electric fields lower than those given by the estimate (j42|. 



III. RESPONSE OF GRAPHENE TO A PULSE EXCITATION 



Now consider the response of graphene to a pulse excitation E'^^*(t) = Eoro(5(i), where Eo and To are the amplitude 
and the duration of the pulse. The pulse amplitude Eq can be arbitrarily strong. Eq. ([M]) now assumes the form 



(43) 

The current is determined again by Eq. ([12]). The solution of Eq. ([43]) can be written in the following implicit form 



Tt = 



-dP, Q = 



2P 



'Pit) PQ{Q) " i + P^' 

where P{t) is the projection of the vector P(t) on the direction of the external electric field and 



Po = P{0) 



PF 



(44) 



(45) 



The function G{Q) is known from Appendix [K\ 
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time, tr time, tr 

FIG. 7: (Color online) The time dependence of (a) the momentum P{t)/Po and (b) the electric current j{t)/ensV, at a pulse 
excitation of graphene. Different curves correspond to different pulse amplitudes Po — cEoto/pf- 

Figure [7^ shows the dependence P{t), calculated from Eq. (H^ . at different values of the electric field parameter 
Po = gEoTq/pf- If the external field is small, Pq ^ 1, the system relaxes after the pulse excitation exponentially, 
similar to the conventional 2D electron systems with the parabolic dispersion, 

P(t) =Poexp(-R), Po<l. (46) 

The characteristic decay time is determined in this case by the inverse radiation decay rate ([35]) . If the external field 
is strong, Pq 3> 1, the momentum of the system first decays linearly in time, 

P{t) = Pa-Tt, Po>l. (47) 

The linear dependence (|47)) remains valid until P{t) reduces down to P{t) ~ 1 (until t ~ Pq/F); after that P{t) decays 
exponentially like in (|46|) . The current j(i) in the strong excitation regime Pq ^ 1 equals its highest possible value 
eUsV and time-independent ai t ^ Po/r? and then exponentially decays (at the time scale ~ F^^) down to zero, as 
shown in Figure [7)3. 

The fact, that after the pulse excitation electrons in graphene move with a constant velocity V ~ 10^ cm/s for quite 
a long time ~ Po/F, may have interesting applications. In a finite-size graphene sample such excited electrons will 
be reflected by the boundaries and oscillate in the sample with the typical frequency V/L, lying in the terahertz 
range, if the sample dimensions L < 1 /im. As at Pq ^ 1 the time Pq/F is much longer than the oscillation period, 
this may lead to a coherent terahertz radiation from graphene excited by a strong pulse electric field. 

IV. SUMMARY AND CONCLUSIONS 

We have developed the quasi-classical kinetic theory of the non-linear electromagnetic response of graphene. Elec- 
trodynamic equations, describing the self-consistent response of the system to a uniform time-dependent external 
electric field, have been derived and solved for the cases of harmonic and pulse excitation. The presented theory is 
valid at hfl ^ max{/i, T}, which covers the frequency range up to ~ 10 — 30 THz for relevant experimental conditions. 

If the system is subjected to a pulse excitation and the amplitude of the external field is small (the linear regime), its 
response is described by a standard exponential decay with the characteristic decay rate F. The density dependence 

1/2 

of the radiative decay rate in graphene (F cx ) differs from that of conventional 2D electron systems {Tpar oc n^). 
In the non-linear regime, when the pulse amplitude is strong, response of the system is not exponential. The average 
momentum of graphene electrons falls down linearly after the pulse excitation, with the average current remaining 
constant, during the time proportional to the amplitude of the external field and inversely proportional to F. 

The graphene layer excited by the harmonic electromagnetic wave with the frequency Q, re-emits radiation at 
higher harmonics mil,, m = 3, 5, 7, . . ., with the efficiency falling down slowly, as l/|m|. The amplitude of the external 
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electric field, required for getting into the non-linear regime, grows with the charge carrier density and is of order of 
several hundred V/cm for typical experimental parameters. The operating frequency of such a frequency multiplier 
can vary in a broad range, from microwaves up to mid-infrared. The efficiency of the frequency upconversion effect 
in graphene can be increased further by using the plasmon, the cyclotron, or the magnetoplasmon resonances. 

The predicted non-linear phenomena in graphene open up new exciting opportunities for building electronic and 
optoelectronic devices for terahertz and sub-terahertz part of the electromagnetic spectrum. 
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APPENDIX A: FUNCTION g{Q) 

The function 

g{Q) ^ f ^ cosOde (^/1 + Qcos0- yJ\-Q cos B] (Al) 
Ti'Q Jo ^ ^ 



can be expanded in powers of Q, 

" 22'-!/!^ + 1)! ■ ^ > 

Since |Q| = |2P/(1 + < 1, the expansion (|A2|) is valid at all values of P. The first terms of this expansion are 

^(g) 1 + (A3) 



If Q = 1, then 



g{\) = ^ « 1.200. (A4) 
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